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Abstract. We consider the continuous and discrete-time Hamilton's variational principle on phase space, 
and characterize the exact discrete Hamiltonian which provides an exact correspondence between discrete 
and continuous Hamiltonian mechanics. The variational characterization of the exact discrete Hamiltonian 
naturally leads to a class of generalized Galerkin Hamiltonian variational integrators, which include the 
symplectic partitioned Runge-Kutta methods. We also characterize the group invariance properties of 
discrete Hamiltonians which lead to a discrete Noether's theorem. 



1. Introduction 



1.1. Discrete Mechanics. Discrete-time analogues of Lagrangian and Hamiltonian mechanics, which are 
derived from discrete variational principles, yield a class of geometric numerical integrators [13] referred to as 
variational integrators [13 [H] . The discrete variational approach to constructing numerical integrators is of 
interest as they automatically yield methods that are symplectic, and by a backward error analysis, exhibit 
(~| . bounded energy errors for exponentially long times (see, for example, [H]). When the discrete Lagrangian 

^ ' or Hamiltonian is group-invariant, they will yield numerical methods that are momentum preserving. 

Discrete Hamiltonian mechanics can be derived from discrete Lagrangian mechanics by relaxing the dis- 
crete second-order curve condition. The dual formulation of this constrained optimization problem yields 
discrete Hamiltonian mechanics jl5|. Alternatively, the second-order curve condition can be imposed using 
Lagrange multipliers, and this corresponds to the discrete Hamilton-Pontryagin principle jl6) . 
. In contrast to the prior literature on discrete Hamiltonian mechanics, which typically start from the 

' Lagrangian setting, we will focus on constructing Hamiltonian variational integrators from the Hamiltonian 

. point of view, without recourse to the Lagrangian formulation. When the Hamiltonian is hyperregular, it 

is possible to obtain the corresponding Lagrangian function, adopt the Galerkin construction of Lagrangian 
variational integrators to obtain a discrete Lagrangian, and then perform a discrete Legendre transformation 
' to obtain a discrete Hamiltonian. This is described in the following diagram: 

H{q,p) — ^L{q,q) 

I 

• ;-( . I 

X 



Hj{qo,Pi) < Ldiqo,qi) 

The goal of this paper is to directly express the discrete Hamiltonian in terms of the continuous Hamiltonian, 
so that the diagram above commutes when the Hamiltonian is hyperregular. An added benefit is that such 
an approach would remain valid even if the Hamiltonian is degenerate, as is the case for point vortices (see 
[22', p. 22), and no corresponding Lagrangian formulation exists. 

The Galerkin construction for Lagrangian variational integrators is attractive, since it provides a general 
framework for constructing a large class of symplectic methods based on suitable choices of finite-dimensional 
approximation spaces, and numerical quadrature formulas. Our approach allows one to apply the Galerkin 
construction of variational integrators to Hamiltonian systems directly, and may potentially generalize to 
variational integrators for multisymplectic Hamiltonian PDEs [H |T71 |TS] . 

Discrete Lagrangian mechanics is expressed in terms of a discrete Lagrangian, which can be viewed as a 
Type 1 generating function of a symplectic map, and discrete Hamiltonian mechanics is naturally expressed 
in terms of discrete Hamiltonians |15j . which are either Type 11 or 111 generating functions. The discrete 
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Hamiltonian perspective allows one to avoid some of the technical difficulties associated with the singularity 
associated with Type I generating functions at time t — (see [19], p. 177). 

Example 1. To illustrate the difficulties associated with degenerate Hamiltonians, consider 

H{q,p) = qp, 

with Legendre transformation given by ¥H : T*Q — > TQ, {q,p) {q,dH/dp) ~ {q,q)- Clearly, in this 
situation, the Legendre transformation is not invertible. Furthermore, the associated Lagrangian is identically 
zero. I.e., L{q,q) = pq - H{q,p)\.^gjjjQ^ ^ pq - qp\q=q = 0. 

The associated Hamilton's equations is given by q — dH/dp — q, p — —dH/dq — —p, with exact solution 
q{t) — q{0)exp(t), p(t) — p(0)exp(— <). This exact solution is, in general, incompatible with the {qo,qi) 
boundary conditions associated with Type I generating functions, but it is compatible with the {qa,pi) boundary 
conditions associated with Type II generating functions. 

In view of this example, our discussion of discrete Hamiltonian mechanics will be expressed directly in 
terms of continuous Hamiltonians and Type II generating functions. 

1.2. Main Results. We provide a characterization of the Type II generating function that generates the 
exact flow of Hamilton's equations, and derive the corresponding Type II Hamilton-Jacobi equation that it 
satisfies. By considering a discrete Type II Hamilton's variational principle in phase space, we derive the 
discrete Hamilton's equations in terms of a discrete Hamiltonian. We provide a variational characterization 
of the exact discrete Hamiltonian that, when substituted into the discrete Hamilton's equations, generates 
samples of the exact continuous solution of Hamilton's equations. Also, we introduce a discrete Type II 
Hamilton-Jacobi equation. 

From the variational characterization of the exact discrete Hamiltonian, we introduce a generalized 
Galerkin approximation from both the Hamiltonian and Lagrangian sides, and show that they are equiv- 
alent when the Hamiltonian is hyperregular. In addition, we provide a systematic means of implementing 
these methods as symplectic-partitioned Runge-Kutta (SPRK) methods. We also establish the invariance 
properties of the discrete Hamiltonian that yield a discrete Noether's theorem. Galerkin discrete Hamilto- 
nians derived from group-invariant interpolatory functions satisfy these invariance properties, and therefore 
preserve momentum. 

1.3. Outline of the Paper. In Section [2l we present the Type II analogues of Hamilton's phase space 
variational principle and the Hamilton-Jacobi equation, and we consider the discrete-time analogues of 
these in Section [S] In Section |4l we develop generalized Galerkin Hamiltonian and Lagrangian variational 
integrators, and consider their implementation as symplectic-partitioned Runge-Kutta methods. In Section 
[51 we establish a discrete Noether's theorem, and provide a discrete Hamiltonian that preserves momentum. 



2.1. Hamilton's Variational Principle for Hamiltonian and Lagrangian Mechanics. Considering 
a n-dimensional configuration manifold Q with associated tangent space TQ and phase space T*Q. We 
introduce generalized coordinates q = {q^,q^, . . . ,q") on Q and {q,p) = {q^,q^, . . . ,q"',pi,p2, . . . ,Pn) on 
T*Q. Given a Hamiltonian H : T*Q — >■ M, Hamilton's phase space variational principle states that 



where extp denotes the extremum over p. Then, Hamilton's phase space principle is equivalent to Hamilton's 
principle. 



2. Variational Formulation of Hamiltonian Mechanics 




for fixed q{0) and q{T). This is equivalent to Hamilton's canonical equations, 
... .OH . dH 

(1) q = -g^[q,p), p^ --Q^[q,p)- 

If the Hamiltonian is hyperregular, there is a corresponding Lagrangian L : TQ R, given by 



L{q,q) = ext pq - H{q,p) = pq - p) l^^a^/ap , 
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for fixed q{0) and q{T). The exact discrete Lagrangian is then given by, 



(90,91) = ext 

qec^{[o.h]..Q) Jo 

q(0)=qo,q{h)=qi 



L{q{t)Ait))dt 



(?,p)6C"([0,/i],T-Q) Jo 

q{o)=qo,q{h)=qi 



PQ - H{q,p)dt, 



which correspond to Jacobi's solution of the Hamilton-Jacobi equation. The usual characterization of the 
exact discrete Lagrangian involves evaluating the action integral on a curve q that satisfies the boundary 
conditions at the endpoints, and the Euler-Lagrange equations in the interior, however, as we will see, the 
variational characterization above naturally leads to the construction of Galcrkin variational integrators. 

2.2. Type II Hamilton's Variational Principle in Phase Space. The boundary conditions associated 
with both Hamilton's and Hamilton's phase space variational principle are naturally related to Type I 
generating functions, since they specify the positions at the initial and final times. We will introduce a 
version of Hamilton's phase space principle for fixed q{0), p{T) boundary conditions, that correspond to 
a Type H generating function, which we refer to as the Type H Hamilton's variational principle in phase 
space. As would be expected, this will give a characterization of the exact discrete Hamiltonian. Taking 
the Legendre transformation of the Jacobi solution of the Hamilton-Jacobi equation leads us to consider the 
following functional, 6 : C'^{[G,T],T*Q) R, 

(2) <S{q{-),p{-))=p{T)q{T)- f [pq - H{q{t),pmdt. 



Lemma 1. Consider the action functional &{q{-),p{-)) given by The condition that &{q{-),p{-)) is 

stationary with respect to boundary conditions Sq{0) — 0, Sp{T) = is equivalent to iq{-),p{-)) satisfying 
Hamilton's canonical equations ([T|). 

Proof. Direct computation of the variation of & over the path space C'^{[0,T],T*Q) yields, 
SG^q{T)Sp{T)+p{T)6qiT) 

q{t)5p{t) +P{t)5q{t) - ^ m.P{t)) 5q{t) - ^ {q{t),p{t)) Sp{t) 

By using integration by parts and the boundary conditions Sq{0) — 0, Sp{T) — 0, we obtain, 
(3) 66 = q{T)5p{T) + p{T)5q{T) - p{T)5q{T) + p{0)5q{Q) 

dH 



dt. 



m 



m 

dH 



{q{t),p{t)) ] Sq{t) - ( q{t) - ^{q{t),p{t))) 5p{t) 



(g(t),p(t)) Ug(t)- 



dH 
dp 



dp 

{q{t)Mt))] 5p{t) 



dt 



dt. 



If (9,^) satisfies Hamilton's equations ([T]), the integrand vanishes, and 5& 
5& = Q for any 5q{Q) = 0, 5p{T) = 0, then from ([3]), we obtain 



0. Conversely, if we assume that 



6& = 



T r 



P{t) + ^m,P{t))) Sq{t) (<ZW - ^iqit),pit)) ] Spit) 



dt = 0, 



and by the fundamental theorem of calculus of variations [2j, we recover Hamilton's equations. 



m 



dH 
dp 



iqit),p{t)) 



m 



dH 
dq 



□ 



The above lemma states that the integral curve {q{-),p{-)) of Hamilton's equations extremizes the action 
functional &{q{-),p{-)) ^ for fixed boundary conditions q{0), p[T). We now introduce the function S{qo,PT), 
which is given by the extremal value of the action functional 6 over the family of curves satisfying the 
boundary conditions g(0) = qo, p(T) = px, 

(4) 5((?o,Pt) = ext ©(q(.),p(.)) = ext prqr - [pq - H{q,p)] dt. 

{q,p)eC^{[0,T].T-Q) (g,p)eC"([0,T],T-Q) Jo 

q(0)=qo.p(T)=PT q(0)=qo,p(T)=PT 
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The next theorem describes how 5((Jo,Pt) generates the flow of Hamihon's equations. 

Theorem 1. Given the function S {qa , pt) defined by the exact time-T flow map of Hamilton's equations 
(QotPo) ^ (QTtPt) is implicitly given by the following relation, 

(5) qr ^ D2S{qo,pT), pa = DiS{qo,PT)- 

In particular, S{qo,pT) is a Type II generating function that generates the exact flow of Hamilton's equations. 
Proof. We directly compute 



dS 
dqo 



{qo,PT) 



dqr 
dqo 

dqr 
dqo 

Po - 



PT 



PT 



^''-^m + '-^Pit) 



dqjt) dH 
dqo dq 



(q^p) 



dqo 
dp{t) 
dqo 



10 



-Po 

10 

dH 
dp 



T r 



il,P) - 



dp{t) 
dqo 
dq{t) 



q 



dqo 



dH 

dp 
dH 
dq 



(q^p) 
iq,p) 



dp{t) dH 
dqo dp 
dq{t) 



{q^p) 



dt 



10 



dH 
dq 



(q^p) 



dt 



dt. 



where we used integration by parts. By Lemma [1] the extremum of © is achieved when the curve {q,p) 
satisfies Hamilton's equations. Consequently, the integrand in the above equation vanishes, giving pq = 
■§^{qo,PT)- Similarly, by using integration by parts, and restricting ourselves to curves {q,p) which satisfy 
Hamilton's equations, we obtain 



dS , . dpT dqr 
— [qo,PT) = t; — qr + — PT 

OpT OpT OpT 



dpit).,, dq(t) dq{t)dH, , dp{t)dH, / 

^ 'qit) + -4^p{t) - -^^{q,p) - -^^{q,p) 



= qr- 
= qT- 



dqj 
dpi 



-PT 



dqi 
dpi 



-PT 



dpT 
dqo 



dp: 



-po 



dpT 
^ ' dp(t) 
dpT 



dpT dq 
dH 



dpT dp 
dpT V 



dt 



dH 



dt 



□ 



2.3. Type II Hamilton— Jacobi Equation. Let us explicitly consider 5(90 7 Pt) as a function of the time T, 
which we denote by 5t(qo:?'t)- Theorem [T] states that the Type H generating function Sxiqo^PT) generates 
the exact time-T flow map of Hamilton's equations, and consequently it has to be related by the Legendre 
transformation to the Jacobi solution of the Hamilton-Jacobi equation, which is the Type I generating 
function for the same flow map. Consequently, we expect that the function SriqoiPT) satisfies a Type H 
analogue of the Hamilton-Jacobi equation, which we derive in the following proposition. 



Proposition 1. Let 

(6) S2iqo,P,t) =Stiqo,p) 



Gxij 

(g,p)GC"([0,t],T*Q) 

g(o)=go, p{t)=p 



p{t)q{t) 



\p{s)q{s) - H{q{s),p{s))]ds 



Then, the function S2{qo,P-,t) satisfies the Type II Hamilton Jacobi equation, 

dS2{qo,P, t) 



(7) 



dt 



op 



Proof. From the definition of S2(qo,P,t), the curve that extremizes the functional connects the fixed initial 
point {qo,Po) with the arbitrary final point {q,p) at time t. Computing the time derivative of S2{qo,Pit) 
yields. 



(8) 

On the other hand, 
(9) 



^ = mq{t) +p{t)m - p{t)m + H{q{t),p(t)). 



dS2 .,^,dS2 dS2 



DISCRETE HAMILTONIAN VARIATIONAL INTEGRATORS 



5 



Equating ^ and (O, and applying ([5]) yields 

(10) -g^=p{tm+H{q{t),p{t)) = H{q{t),p{t)) ^ H [-q^.p) ■ 

□ 

The Type II Hamilton-Jacobi equation also appears on p. 201 of [T3^ and in [8]. However, this equation 
has generally been used in the construction of symplectic integrators based on Type II generating functions 
by considering a series expansion of 52 in powers of t, substituting the series in the Type II Hamilton-Jacobi 
equation, and truncating. Then, a tcrm-by-term comparison allows one to determine the coefficients in the 
series expansion of S2 , from which one constructs a symplectic map that approximates the exact flow map 

However, approximating Jacobi's solution on the Lagrangian side, or the exact discrete right Hamiltonian 
S{qt)TPT) in Oj in terms of their variational characterization provides an elegant method for construct- 
ing symplectic integrators. In particular, this naturally leads to the generalized Galerkin framework for 
constructing discrete Lagrangians and discrete Hamiltonians, which we will explore in the rest of the paper. 

In Section [3l we will also present a discrete analogue of the Type II Hamilton-Jacobi equation, which 
can be viewed as a composition theorem that expresses the discrete Hamiltonian for a given time interval 
in terms of discrete Hamiltonians for the subintervals. This can be viewed as the Type 11 analogue of the 
discrete Hamilton-Jacobi equation that was introduced in [lO] . 

3. Discrete Variational Hamiltonian Mechanics 

3.1. Discrete Type II Hamilton's Variational Principle in Phase Space. The Lagrangian formula- 
tion of discrete variational mechanics is based on a discretization of Hamilton's principle, and a comprehen- 
sive review of this approach is given in [21]. The Hamiltonian analogue of discrete variational mechanics 
was introduced in 1151, wherein discrete Lagrangian mechanics was viewed as the primal formulation of a 
constrained discrete optimization problem, where the constraints are given by the discrete analogue of the 
second-order curve condition, and dual formulation of this yields discrete Hamiltonian variational mechan- 
ics. An analogous approach is based on the discrete Hamilton-Pontryagin variational principle |16j , in which 
the discrete Hamilton's principle is augmented with a Lagrange multiplier term that enforces the discrete 
second-order curve condition. 

We begin by introducing a partition of the time interval [0, T] with the discrete times Q = < ti < . . . < 
= T, and a discrete curve in T*Q, denoted by {{qk,Pk)}^^Q, where qk ~ q{tk)^ and pk ~ p{tk)- Our 
discrete variational principle will be formulated in terms of a discrete Hamiltonian {qk,Pk+i): which is 
an approximation of the Type II generating function given in 

ntk+i 

(11) H'^iqk^Pk+i) ~ ext p{tk+i)q{tk+i) - j [pq - H{q,p)]dt. 

iq,p)eC^iltk,tk+l],T'Q) Jtk 
qitk)=qk,p{tk + l)=Pk + l 

As we saw in Section [2l the curve in phase space with fixed boundary conditions (qo^px) that extremizes 
the functional ([2]), 

&{q{-):P{-)) = p{T)q{T) ~ Clpq - H{q{i),p{t))]dt, 

Jo 

satisfies Hamilton's canonical equations. Consequently, we can formulate discrete variational Hamiltonian 
mechanics in terms of a discrete analogue of this functional, which is given by, 

^-1 rtk+i 

(12) ed{{{qk,Pk)}Lo) = PNqN - E / N - H{q{t),p{t))] dt 

k=Q 

= PNqN - ^ [pk+iqk+1 - Hj{qk,Pk+i)] ■ 

fc=0 

Then, the Type II discrete Hamilton's phase space variational principle states that ^©d({('Zfe,Pfe)}fcLo) ~ 
for discrete curves in T*Q with fixed {qo,pN) boundary conditions. 
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Lemma 2. The Type II discrete Hamilton's phase space variational principle is equivalent to the discrete 
right Hamilton's equations 

qk = D2H+iqk-i,Pk), k = l,...,N-l, 

\ / -I- 

Pk = DiHJ{qk,pk+i), k = l,...,N -1, 

where {qk,Pk+i) is defined in pT|) . 

Proof. We compute the variation of (3d, 

/ JV-l \ 

5&d = 5 \pNqN - ^{Pk+iqk+i - Hj {qk,Pk+i)) 1 

(N-2 N~l \ 

- X! Pk+iqk+i + ^ H+{qk,Pk+i) I 
k=0 fc=0 / 

= - X! i<lk+iSpk+i +Pk+iSqk+i) + X! {DiH^{qk,Pk+i)5qk + D2Hj {qk,Pk+i)Spk+i) 

k=0 k=0 
N-1 N-1 

= - '^{qkSpk + PkSqk) + ^ DiHj{qk,Pk+i)5qk + DiHj{qQ,pi)Sqo 

k=l k=l 
N-1 

+ ^ D2H^{qk-i,Pk)Spk + D2H^{qN~i,PN)SpN 

k=l 

N-1 N-1 

= - ^ {(Ik- D2Hj{qk-i,Pk)) Spk {Pk- DiHj{qk,Pk+i)) Sqk 

k=l k=l 

+ DiH^{qo,pi)dqo + D2H^{qN-i,PN)SpN- 

where we reindexed the sum, which is the discrete analogue of integration by parts. Using the fact that 
(qotPn) are fixed, which imphes Sqo — 0, Sp^ = 0, the above equation reduces to 

N-1 N-1 

(14) <56d = - ^ {qk - D2H+{qk-i,Pk)) Spk ~Y.{P^- DiH+{qk,Pk+l)) Sqk. 

k=l k=l 

Clearly, if the discrete right Hamilton's equations, qk = D2H'^ {qk-i,pk), Pk = DiHd ilk,Pk+i), are satisfied, 
then the functional is stationary. Conversely, if the functional is stationary, a discrete analogue of the 
fundamental theorem of the calculus of variations yields the discrete right Hamilton's equations. □ 

The above lemma states that the discrete-time solution trajectory of the discrete right Hamilton's equa- 
tions (fT3)) extremizes the discrete functional (fT2|) for fixed qQ,PN- However, it does not indicate how the dis- 
crete solution is related to po, q^- Note that the discrete solution trajectory that renders &d{{{{qk,Pk)}k=o) 
stationary depends on the boundary conditions go, Pn- Consequently, we can introduce the function Sd 
which is given by the extremal value of the discrete functional Gd as a function of the boundary conditions 
q{to),p{tj\j), and is explicitly given by 

(15) Sd{q{to),pitN))= ^ ext ^ 6d{{{qk,Pk)}Lo) 

qo=q(to),PN=p(tN) 

N-1 

= ext pNqN -^{Pk+iqk+l- H+{qk,Pk+l))■ 

{qk■,Pk)^T■'^Q ^ 

9o=<3(*o),P]v=p(tjv) 

Then, using a similar approach to the proof of Theorem [l] we compute the derivatives of Sd{qo,PN) with 
respect to qo,pN- By reindexing the sum, which is the discrete analogue of integration by parts, we obtain 

N-2 N-1 \ 



^i<10,PN) = ^ ( - 51 Pk+lQk+l + Hj{qk,Pk+l) 
" \ k=0 k=0 / 
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= - V ^ (<Zfc - D2H+iqk^,,pk)) - V Tp^ (Pfc - D,H+{qk,Pk+i)) + D,H+iqo,pi). 

By Lemma [51 the extremum of &d is obtained if the discrete curve satisfies the discrete right Hamihon's 
equations (|13p . Thus, by the definition of Sdiqa^PN), the above equation reduces to 

(16) DiSd{qQ,PN) = DiHj{qQ,pi). 

A similar argument yields 

(17) 

dSd 

dpN 



N-2 

{qQ,PN) = I - V Pk+iQk+i 





= ~ X! ^ D2H^{qk-i,Pk)) - X! (^^^^ ^ DiHj{qk,Pk+i)) + D2Hj {qN-i,PN) 

fc— 1 fc— 1 

= D2Hj{qN^l,PN)- 

Recall that the exact discrete Hamiltonian S{qo,pT) defined in Q is a Type II generating function of the 
symplectic map, implicitly defined by the relation ([5]), that is the exact flow map of the continuous Hamilton's 
equations. To be consistent with this, we require Sd{qo,PN) satisfies the relation ([5]), which is to say 

(18) qN ^ D2Sd{qo,PN), Po = DiSd{qo,PN)- 
Comparing (HHl-dTTl) and (HH]) we obtain 

(19) qN = D2Hj{qN-i,PN), po ^ DiHj {qo,pi). 

Then, by combining (|13p and (IT9| . we obtain the complete set of discrete right Hamilton's equations 

(20a) qk+i=D2H+{qk,Pk+i), fc = 0, 1, . . . , iV - 1, 

(20b) pk = DiH+{qk,Pk+i), fc = 0,l,...,iV-l. 

It is easy to see that 

= ddHj{qk,Pk+i) = d(DiHj{qk,Pk+i)dqk + D2Hj {qk,Pk+i)dpk+i) 

= d {pkdqk + qk+idpk+i) = dpk A dqk - dpk+i A dqk+i, 

for k = 0, . . . , N — I. Then, successively applying above equation gives 

dpa A dqa = dpi A dqi = ■ • • = dpN-i A dq^-i = dpM A dqN ■ 

This implies that the map from the initial state (qojPo) to the final state {qN,PN) defined by (jl8p is sympletic, 
since it is the composition of N symplectic maps {qk,Pk) ^ {ik+iiPk+i), k = 0, . . . , N ~ 1, which are given 
by ([20a| -(|20b |) . Ahernatively, one can directly prove symplecticity of the map {qo,Pa) (qntPn) by using 
(fTSl) to compute = d'^Sd{qo,PN) — dpQ A dqa — dpN A dqN- Given initial conditions qo^Po, and under the 

regularity assumption "q^^q^^ -{qk,Pk+i)\ 7^ 0, we can solve (|20b[) to obtain pi, then substitute pi into pOa[) 

to get qi. By repeatedly applying this process, we obtain the discrete solution trajectory {{qk,Pk)}k=i- 

3.2. Discrete Type II Hamilton— Jacobi Equation. A discrete analogue of the Hamilton- Jacobi equa- 
tion was first introduced in [10], and the connections to discrete Hamiltonian mechanics, and discrete optimal 
control theory were explored in [23' . In essence, the discrete Hamilton- Jacobi equation therein can be viewed 
as a composition theorem that relates the discrete Hamiltonians that generate the maps over subintervals, 
with the discrete Lagrangian that generates the map over the entire time interval. 

We will adopt the derivation of the discrete Hamilton-Jacobi equation in [23', which is based on intro- 
ducing a discrete analogue of Jacobi's solution, to the setting of Type II generating functions. 



Theorem 2. Consider the discrete extremum function (|15|) ." 

fc-i 

(21) <Sd{Pk) = PkQk - V [pi+iqi+i - H+{qi,pi+i)] 



1=0 
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which can be obtained from the discrete functional (|12p by evaluating it along a solution of the right discrete 
Hamilton's equations Ij20p . Each Sll{pk) is viewed as a function of the momentum pk at the discrete end 
time ifc. Then, these satisfy the discrete Type II Hamilton-Jacobi equation: 

(22) - Si{pu) - H+{DS^^{pu).pu+i) - Pk ■ DS^ipk). 
Proof. From Eq. (|21|) . we have 

(23) S''+\pk+i) - Si(pk) = H+{qk,Pk+i) - Pk ■ gfe, 

where qk is considered to be a function of pk and pk+i, i-c, qk — qk{pk,Pk+i)- Taking the derivative of both 
sides with respect to pk , we have 

-DS'^M = -qk + ^ ■ [DiH+{qk,pk+i)-Pk] . 

Opk 

However, the term in the brackets vanish because the right discrete Hamihon's equations ((^(7)) are assumed 
to be satisfied. Thus we have 

(24) qk = DS^ipk). 

Substituting this into ^ gives ([22]). □ 



3.3. Summary of Discrete and Continuous Results. We have introduced the continuous and discrete 
variational formulations of Hamiltonian mechanics in a parallel fashion, and the correspondence between the 
two are summarized in Figure [TJ Similarly, the correspondence between the continuous and discrete Type II 
Hamilton-Jacobi equations are summarized in Table [T] 



Figure 1. Continuous and discrete Type II Hamilton's phase space variational principle. 
In the continuous case, the variation of the action functional © over the space of curves 
gives Hamilton's equations, and the derivatives of the extremum function <S with respect to 
the boundary points yield the exact flow map of Hamilton's equation. In the discrete case, 
the variation of the discrete action functional &d over the space of discrete curves gives 
the discrete right Hamilton's equations and the derivatives of extremum functional Sd with 
respect to the boundary points yield the symplectic map from the initial state to the final 
state. 
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pQ ^ _Di5(go,pT) 



Disc. Action 
Functional 



Disc. Right 
Hamilton's Eqn. 

pk ^ DiH'^ {qk,pk+i) 



Disc. Extremum 
Function Sd 



Qn ^ D2Sd{qo.pN) 
po ^ DiSd{qo,PN) 



Symplccticity 
— ddS — dpo A 

dqo — dpT A dqT 



Symplccticity 
— ddSd — dpo A 
dqo — dpN A dqN 
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Table 1. Correspondence between ingredients in the continuous and discrete Type II 
Hamilton-Jacobi theories; No is the set of non-negative integers and ]R>o is the set of 
non-negative real numbers. 



Continuous 


Discrete 




\Pk, Kj fc (i^ X i"*lo 


q — dH/ dp, 
p = —dH/dq 


qk+1 — ^2J^d (qk,Pk+i), 
Pk = DiH+{qk,Pk+i) 


S2{p, t) = p(t)q{t) - f [p{s)q{s) - H{q{s),p{s))] ds 
Jo 


fc-i 

SdiPk) = Pkqk - ^ [pi+iqi+i - Hj{qi,pi+i)] 

1=0 


db2 = dp + — — dt 
dp at 


S^+\p,+,) - S^Pk) 


qdp + H{q,p) dt 


Hai<lk,Pk+i) -Pk-qk 


dS2 rr(dS2 \ 


'Sd^^(Pfc+i) -Sa{pk) 

= H+(DSl{pk),Pk+i)-Pk ■ DS^ipk) 



4. Galerkin Hamiltonian Variational Integrators 

4.1. Exact Discrete Hamiltonian. The exact discrete Hamiltonian ^^j~exact(90 7Pi) and the discrete right 
Hamilton's equations ()20p generate a discrete solution curve {qk,Pk}^=o that samples the exact solution 
{q{-),p{-)) of the continuous Hamilton's equations for the continuous Hamiltonian H(q,p), i.e., qk = q{tk) 
and Pk = p{tk)- 

By comparing the definition of a discrete right Hamiltonian function H^{qQ,pi) on [0,h] and the 
corresponding discrete Hamiltonian flow in (|20p to the definition (|4|) of the extremal function on [0,T] and 
corresponding symplectic map given by ([5]) , and applying Theorem [TJ it is clear that the discrete right 
Hamiltonian function on [0, /i], given by 

(25) H+ iqo,p,)^ ext pm- [p{t)m-H{q{t),pmdt, 

(q.p)eC^{[0,T],T'Q} Jq 

q{to)=qo,p{ti)=pi 
is an exact discrete right Hamiltonian function on [0,h]. 

4.2. Galerkin Discrete Hamiltonian. In general, the exact discrete Hamiltonian is not computable, 
since it requires one to evaluate the functional (3(g(-),p(-)) given in ([2]) on a solution curve of Hamilton's 
equations that satisfies the given boundary conditions {qo,pi). However, the variational characterization of 
the exact discrete Hamiltonian naturally leads to computable approximations based on Galerkin techniques. 
In practice, one replaces the path space C^([0, T],T*Q), which is an infinite-dimensional function space, with 
a finite-dimensional function space, and uses numerical quadrature to approximate the integral. 

Let {^/'i(r)}f^i, T G [0, 1], be a set of basis functions for a s-dimensional function space Cj. We also choose 
a numerical quadrature formula with quadrature weights bi, and quadrature points Ci, i.e., f{x)dx ~ 
From these basis functions and the numerical quadrature formula, we will systematically 
construct a generalized Galerkin Hamiltonian variational integrator in the following manner: 

1. Use the basis functions ^pi to approximate the velocity q over the interval [0,h], 

s 

(jd(T/i) =^FV,;(t). 
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2. Integrate qd{p) over [0, r/i], to obtain the approximation for the position g, 



qdirh) 



Jo Jo 



where we apphed the boundary condition qd{0) ~ qo- Applying the boundary condition qd{h) = qi 
yields 

qi = qd{h) - go + V / Mp)dp ^ qo + hY,B^V\ 

z=l ^=l 

where Bi = ijji(T)dT. Furthermore, we introduce the internal stages, 

Q' = qd(cih) ^qo + hy^V^ / ^3 {T)dT = go + /^ V AijV , 
where Aij — J^^ ^pj{T)dT. 

3. Let = p{cih). Use the numerical quadrature formula {bi,Ci) and the finite-dimensional function 
space Cj to construct H^{qo,pi) as follows 

.h 

H+{qo,Pi)- ext piqi- [p{t)q{t)-H{q{t),p{t))]dt 



{q,p)ec^[[o,T\.T'Q) 

q(to)=qo,p{tl)=Pl 



(26) 



Hd{qo,Pi)= ext lpiqd{h)-hy^bi[p{cih)qd{cih)-H{qd{cih),p{cih))] 

gd&C'([o,h]'Q),Pi<^Q'' ^ 



^ext _ |pi l^qo + ^ j - ^ h, 

ext K{qo,V\P\pi). 



To obtain an expression for H'^ {qo , pi ) , we first compute the stationarity conditions for K{qQ , X^* , P* , pi ) 
under the fixed boundary condition {qo,pi). 

(27a) = ^-^^(gO'^^^^/'^'Pi) ^ hp^B, -hj^h [P'^c,) hA,~{Q\n^ , j = 

j = l,...,s. 



z— 1 ^ 

(27b) . (i:*(«-)v" - f (Q^^-')) . 

4. By solving the 2s stationarity conditions (P7|. we can express the parameters V^,P^, in terms of 
go, Pi, i-e., 1^* = y (go,Pi) and = P*(qo,Pi)- Then, the symplectic map (<7o,Po) H> (<?i,pi) can be 
expressed in terms of the internal stages 

dH+{qo,pi) _ dK{qo,V''{qo,pi),P'{qo,pi),pi) 



(28) po 



dqo dqo 
dK dK dV dK dP' dK 



dqo dV' dqo dP' dqo 

Pi + hy2b,—iQ\p^), 



Similarly, we obtain 
(29) qi 



dH+{qo,p,) ^ dK{qo,VHqo,Pl),P'{qo,Pl),Pl) 

dpi dpi 
dK dV dK dP' dK _ dK 
dV- dpi dP^ dpi dpi dpi 
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^qo + hJ2B^V\ 

1=1 

Without loss of generality, we assume that the quadrature weights 6, ^ 0. Then, the stationarity condition 
(|27bp reduces to J2Ui i^^i^jW' - ^{Q^,P^) = 0. Moreover, by substituting ^ into the stationarity 
condition (HZil), we obtain, Y.Ui hP'^j{ci) - p^B, + hJ2U^{b^Bj - h,A,j)^{Q\P') = 0. 

In summary, the above procedure gives a systematic way to construct a generalized Galerkin Hamiltonian 
variational integrator, which can be rewritten in the following compact form, 

(30a) qi=qa + hY,B,V\ 

1=1 

(30b) p,=p,-hY,h-g-{Q\P'), 

i—1 



(30c) Q' = qo + hJ2^^Jy', 1 = 1,.. 



(30d) = Y,b^P'^Ac^)~PoB,+hJ2{b^BJ~hA,,) — {Q\P'), j = l,.. 



1=1 



(30e) o = Y,Mcj)V' --g-iQ',P'), J = 

1=1 ^ 

where (&i, q) are the quadrature weights and quadrature points, and Bi = ^pi{T)dT, Aij = J^' ^pj(T)dT. 

This is the general form of a Galerkin Hamiltonian variational integrator. Issues of solvability, convergence, 
and accuracy, depend on the specific Hamiltonian system, and the choice of finite-dimensional function space 
and numerical quadrature formula {bi,Ci). We will not perform an in depth analysis here, but we will 
illustrate how our proposed framework is related to the discrete Lagrangian based methods given in [21] and 
p. 209 of 

4.3. Galerkin Variational Integrators from the Lagrangian Point of View. In this subsection, 
we investigate the generalized Galerkin variational integrators from Lagrangian point of view when the 
Hamiltonian function is hyperregular. In this case, the exact discrete right Hamiltonian function is related 
by the Legendre transformation to exact discrete Lagrangian function, i.e., 

(31) Lr'\qo,qi)= ext f L{q{t)A{t))dt 

9(0)=9o.g('0='?i 



ext / pq — H{q,p)dt 



iq,p)<£C'ao,h],T'Q) 
qiO)=qo,q{h)=qi 

Pm- H+ ^{qo,pi) 



qi=D2H+ ^{qo,pi) 



We wish to see how Galerkin variational integrators that are derived from the Hamiltonian and Lagrangian 
sides are related. In order for the comparison to make sense, we will approximate the exact discrete La- 
grangian using the same basis functions and numerical quadrature formula as on the Hamiltonian side. As 
before, let {V'i(''')}f=ii ^ [Oj1]j be a set of basis functions for a s-dimensional function space CJ, and 
choose a numerical quadrature formula with quadrature weights bi, and quadrature points c^. From these 
basis functions and the numerical quadrature formula, we will systematically construct a generalized Galerkin 
Lagrangian variational integrator in the following manner: 

1. Use the basis functions ipi to approximate the velocity q over the interval [0,h], 



(r/.)=^FV,:(r). 
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2. Integrate qd{p) over [0, r/i], to obtain the approximation for the position g, 

qd{rh)=qM+ y2^'Mp)d{ph) = qo + hTV' r Mp)dp, 
Jo Jo 

where we applied the boundary condition gd(0) = qo- In the discrete Lagrangian framework, the 
boundary conditions are given by {qo, qi), so we will use a Lagrange multiplier to enforce the boundary 
condition qdih) — qi, 

s „i s 

qi = qd{h) =qa + h^V' / Mp)dp = qo + h^B,V\ 
where Bi — ipi{T)dT. Furthermore, we introduce the internal stages, 
(32) Q'' = qd{c,h) = qo + hY,V' ^j{r)dT = qo + /i^ 



where Aij — J^^ ^/jj{T)dT, and their velocities. 



(33) = gd(c,/^) -^^,(c,)F^'. 

3. Use the numerical quadrature formula (bi,Ci), and the finite-dimensional function space CJ to con- 
struct Ld(qo,qi) as follows 



Ld{qo-,qi) 
(34) Ld{qo,qi) 



ext 

qec^ ilo, h],Q),\ 
■z(o)=go 



cxtj 

9deC^([0,/i],Q),A 
«(o)=go 



L{q{t),q{t))dt 



+ \{qi-q{h)) 



h^b.iL{qd{c.ih),qd{cth)) 



i=l 



ext 

V'.X 



hY^b,L\ qo + hY,A^,V',Yl ^'V-.Xq) 



A(gi - qd{h)) 

^\(qi-qo-hY^B,V' 



i=l 



EE ext K{qo,V\X,qi). 

To obtain an expression for Ld{qo, qi), we first compute the stationarity conditions for K{qo, V^, A, qi) 
under the fixed boundary condition {qo,qi)- 

(35a) = 9K{qo,V\X,q^) ^ j^j^h {^{Q\ Q')hAj + ^{Q\ Q')^3{^^)] ' h\B,, j = 1, . . . , s. 



(35b) = 



dK{qo,V\X,qi) 

dX 



dq 



^qi-qo-hY^B,V' 



4. By solving the 2s stationarity equations (|35|) . we can express the parameters y , A, in terms of qo,qi, 
i.e., = V^{qo,qi), X = X{qo,qi). Then, the symplectic map {qo,Po) ili^Pi) can be expressed in 
terms of the internal stages and the Lagrange multiplier 



(36) 



Po 



dLd{qo,qi) ( dK{qo,V\qo,pi),X{qo,pi),qi) 



dqo \ 
dK dK dV dK dX 



'io 



dqo dV^ dqo dX dqo 

s 

hY,hLq{Q\Q') + X. 



dK 
dqo 



i=l 
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Similarly, we obtain 

dLd{qQ,qi) _ dK{qa,V''{qo,pi),X{qo,pi),qi) 



(37) pi 



dqi dqi 
dK dK dV dK d\ dK 



dqi dV^ dqi dX dqi dqi 
^ A. 

By combining ([5^ and ([57]) . we obtain pi = Pq + hJ2i=i biLq{Q^ ,Q^). Substituting this into the station- 
arity condition yields J^Li hdL / dq{Q\Q^)i^j{ci) - p^Bj - hJ^Uii^Bj - b,A,j)dL/dqiQ\ Q') = 0. 

In summary, the above procedure gives a systematic way to construct a generalized Galerkin Lagrangian 
variational integrator, which can be written in the following compact form, 

s 

(38a) ql^qo + hY,B^V\ 

(38b) Pl^Po + hY,b^-g-(,Q\Q'), 

i—1 

s 

(38c) Q' = qo + hY,A^JV', i = l,...,s 

^ dr dT 
(38d) = hj7{Q\ mA^^) ~ PoBj - hY^ihB, - b,A,) — iQ\Q^), j = 1, . . . , s 

i—1 i—1 

s 

(38e) = Y,Mcj)V' -Q', J = l,.-.,s 

4=1 

where (bi, Ci) are the quadrature weights and quadrature points, and Bi = ipi{T)dT, Ay — J^' il)j{T)dT. 

As expected, this is equivalent to the generalized Galerkin Hamiltonian variational integrator, as the 
following proposition indicates. 

Proposition 2. // the continuous Hamiltonian H{q,p) is hyperregular, and we construct a Lagrangian 
L{q,q) by the Legendre transformation, then the generalized Galerkin Hamiltonian variational integrator 
pOal) - (j30ep and the generalized Galerkin Lagrangian variational integrator (I38ap - (j38ep . associated with the 
same choice of basis functions and numerical quadrature formula {bi,Ci), are equivalent. 

Proof. Since we chose the same basis functions and numerical quadrature formula for both methods, the 
approximations for qi and are the same in both methods, as can be seen by comparing (|30a|) and p8al) . 
pOcp and ([35c]). Since we assumed that the Lagrangian and Hamiltonian are related by L(q, q) = pq—H{q,p), 
subject to the Legendre transformation q = dH/dp{q, p), we consider p to be a function of (q, g), and compute 
dL, .dp dH , , dH , dp dH , , 

-d-q^^^""^ -^^-d-q - 9^^*^'^) ~ -d^^'^^^'^d-q = " 9^^'^'^^' 

dL . . dp dH dp 

^(?,9) = 9-^+p-^(<Z,2>)^=p. 

Since these identities have to hold at the internal stages, we have that 

n TT 

dT 

for i — l,...,s. Clearly, substituting these identities into (j30bp . (j30dp . and (I30ep . yields (j38bp . (j38dP . 
and p8ep . As such the two systems of equations, P0al) - (|30ep and (|38ap - (|38ep . are equivalent, once the 
Legendre transformation and the identities relating the continuous Lagrangian and Hamiltonian are taken 
into account. □ 



14 



MELVIN LEOK AND JINGJING ZHANG 



4.4. Variational integrators and Symplectic Partitioned Runge-Kutta methods. In this subsec- 
tion, we consider a special class of Galcrkin variational integrators, and demonstrate that they can be 
implemented as symplectic partitioned Runge-Kutta methods. 

Let CJ([0, 1], Q) be a s-dimensional function space, and consider a set of basis functions ipiir) on [0, 1], 
and a set of control points c^, i = 1, . . . , s. We would like to construct a new set of basis functions ^^(t) that 
span the same function space, and satisfies 0i(cj) = dij^ where Sij is the Kronecker delta. This is possible 
whenever the matrix 



(39) 



M = 



-01(01) -01(02) 

-02(ci) -02(C2) 



■01 (Cs)' 
-0s (Cs). 



is invertible. In particular, let 0(-) = [-0i(-), . . . ,0s(')]"^i ^-nd construct a new set of basis functions 0(-) 



l(-), . . . , 0s(-)] by 0(-) = M V(')- It is easy to see that 0i(cj) = (5; 



since 



(40) 



01 (Cl) 01 (C2) 
02(ci) 02(C2) 



Cl 



0s (C2) 



01 (Cs) 
02 (Cs) 

0s (Cs) 



1 
1 





We can construct a numerical quadrature formula that is exact on the span of the basis functions ipi (r) as 
follows: Since (j^iicj) = Sij, we can interpolate any function /(r) on [0, 1] at the control points Ci by taking 
/(r) = X]i=i /(ci)0i(''')- Then, we obtain the following quadrature formula. 



(41) 



f{r)d7 



f{T)dT ^ J2 /(c.)0.(T)dT = J2 /(^^) 



1=1 



where bt = (j)i{T)dT are the quadrature weights. By construction, the above quadrature formula is exact for 
any function in the s-dimensional function space CJ([0, 1], Q). Now, if we apply this quadrature formula with 
quadrature points Ci, we obtain a Galerkin variational integrator that can be implemented as a symplectic 
partitioned Runge-Kutta (SPRK) method. 

Theorem 3. Given any set of basis functions ?/'(•) = ['0i('); ■ • ■ j '0s(')]"^; ^^^^ span C^([0,1],Q), consider 
the quadrature formula given in (|41[) . Then, the associated generalized Galerkin Hamiltonian variational 
integrator (j30ap - (|30e[) . which is expressed in terms of the discrete right Hamiltonian function (|26p . can 
be implemented by the following s-stage symplectic partitioned Runge-Kutta method applied to Hamilton's 
equations ([T]); 



(42a) 
(42b) 
(42c) 
(42d) 



qi^qo 



Pi =Po 



Q' = qo 



i=l 



dp 



p^ 



PO 



h^~a,~{Q\P^), 



dq 



.,s. 



.,s. 



where hi = (j)i(T)dT ^ 0, fljj = J^^' <f>j{T)dT, and a. 



bibi 



The basis functions satisfy <j)i (cj ) 



"^3 ~ bi 

and are given by 0(-) = M~^%1]{-), where 0(-) — [0i(-)i ■ • ■ i 0s(-)]'^) '^"'^ defined in ((39l 



Proof. The new basis functions 0(t) are constructed from the original basis functions 0'(t) by the relationship 
0(-) = M~^'il){-). Thus, we have that 0(-) = M(j){-), and in particular, 0'i(r) = X)j=i 0'i(cj)0j(-'')- By 
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substituting this into Bt = 'iJji{T)dT, equation (|30a|) becomes 

s 

qi = qo + hY^B,V' 



1=1 

pi s 



S .IS 

= qo + hJ2v' I] (Cj )0j (r)dT 

^=l j = l 

= qo + hy2—{Q^pn / ^jiT)dT 

i=i ^ 

where we used equation pOeP to go from the third equality to the fourth one. Similarly, by substituting 
'fpkiT) = ^j^i '^k{cj)4>j{T) into Aik = Jq' tl}k{'T)dT and using equation pOep . equation (I30cp becomes 

s 

k=l 

= qQ + h^V^ iJk{cj)4>j{T)dT 



pi dp Jo 

7 = 1 



where a^- = (j)j{T)dT. Note that equation (j30bp has the same form as ()42bp . so we only have to re- 
cover equation (l42dl) . Let E'l' = [sf , . . . , EfY' , where = Ei=i ^^-P'V'ilci) - PoBj + hJ^UiibiB, - 
biAij)^{Q^, F*) 0, which corresponds to pOdp . By substituting Bi — tpi{T)dT — J2'j=i '4'i{cj)4'j{T)dT, 
and bi = (j)i{T)dT into , we obtain 

= ^/ = E - Po^j + /i Y.^hBj - h,A,,) — {Q\ P') 

'i=i i=i 

,= 1 ^=l 

+ VVj(cfc)0fc(r)dr-&, r y2i^jick)Mr)dT] ^{Q\P') 

,:=i V ^" fc=i fc=i / '^'^ 

We swapped the role of the indices i and k in the second to last line to obtain the final equality. Let 
E-^ = [Ef, . . . , where Ef = b,P' - b,po + hJ2l=iihbk - buakr)^{Q\ P""). Then, the above equation 

can be viewed as the j-th component of the system of equations ME'^ = E^ = 0, where M = [V'i(cj)] is 
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invcrtible. Therefore, we have that E'f> = 0, i.e., Ef = 6,F* - hpa + ^ELi(^»^fc " bkaki)^{Q'" , = 



Since 6^ 7^ 0, dividing by bi and recalhng that a,j — 



bibi 



yields (|42d|. 



□ 



Comparison with Discrete Lagrangian SPRK Methods. Proposition [2] states that for hyperregular 
Hamihonians, if one chooses the same basis functions and quadrature formula, the generalized Galerkin 
Haniiltonian variational integrator is equivalent to the generalized Galerkin Lagrangian variational integrator. 
Therefore, the above theorem also applies in the Lagrangian setting. In particular, if one chooses the 
Lagrange polynomials associated with the quadrature nodes Ci as our choice of basis functions i/jiij'), then 
the coefficients of the SPRK method derived above agree with the method derived in [3T] using discrete 
Lagrangians. However, our approach remains valid in the case of degenerate Hamiltonians, for which it is 
impossible to obtain a Lagrangian and apply the method in |21| to derive Hamiltonian variational integrators. 

The derivation on p. 209 of the book [13!, which is analogous to the result in |26], generalizes the approach 
in [21] by considering discrete Lagrangian SPRK methods without the restriction that the Runge-Kutta 
coefficients are obtained from integrals of Lagrange polynomials. It is however unclear how one should 
choose these coefficients. In contrast, our approach provides a systematic means of deriving the coefficients 
by an appropriate choice of basis functions and quadrature formula. Our discrete Hamiltonian method is 
expressed in terms of Type II generating functions and the continuous Hamiltonian as opposed to the discrete 
Lagrangian approach based on Type I generating functions and the continuous Lagrangian. 



Discrete Hamiltonian associated with Galerkin SPRK Method. For the symplectic partitioned 
Rungc-Kutta method (|42p described above, we can explicitly compute the corresponding Type II generating 



function H^{qo,pi) given in ([^ as follows. 



(43) 



H+{qo,pi)=piqdih) -hY^h [P'qMh)-H{Q\P')\ 

s 



= Pi[qo + hJ2b^^{Q\P') 



dH , 
dp 



r) M 

P' — {Q\P')-H{Q\P') 



dH 



piqo 



piqo 



h J2 b^iPi - P1-Q-{Q\P') + hJ2 WQ\ P') 

1=1 P i=l 

^ r)H r)H ^ 

E hci,, — {Q\P')^{Q\P=) + hY,hH{Q\P^). 



dp 



This Type II generating function is consistent with the Type I generating function for SPRK methods that 
was given in Theorem 5.4 on p. 198 of [13]. 

Sufficient Condition for Consistency of the SPRK Method. If the constant function f{x) = 1 is in the 
finite-dimensional function space C^, then by interpolation, we have that 1 = X]i=i f {ci) 4' %{''') — Si=i 
Thus, X]i=i bi — Jq X]i=i 'Pi{'^)dT = 1. Partitioned Runge-Kutta order theory [5] states that the condition 
Xli=i bi = \ implies that the variational integrator ([^^ is at least first-order. Therefore, to obtain a consis- 
tent method, it is sufficient that the constant function is in the span of the basis functions we choose. In 
particular, if we let iPi{t) = 1, we ensure that our method is at least first-order. 

Construction of the SPRK Tableau. Let the symplectic partitioned Runge-Kutta method (|42l) be de- 
noted by the tableau. 



Cl 


an • 


■ ais 


Cl 


an ■ 


■ ais 


Cs 


ftsi • 


(^ss 


Cs 


CLsl ■ 


ass 




bi ■ 


■ bs 




bi ■ 


■ bs 



Based on the above generalized Galerkin method, the coefficients in the partitioned Runge-Kutta tableau 
can be constructed in the following systematic way: 
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. , s. Ensure that M = [ijji{cj)] is invcrtiblc. 

.bs]^ contain the coefRcients in the SPRK tableau. 



Generalized Galerkin Hamiltonian SPRK Method. 

1. Choose a basis set 'tpi{T),T G [0, 1], i = 1, . . . , s, with ■!/'i(r) = 1. 

2. Choose quadrature points C;,?' = 1, 

3. Let the column vector 6 = &2, • ■ 
There are two ways to obtain b: 

i. Compute Bi = %pi{T)dT and let B = [Bi, B2, . . . , B^f. Then, b = M'^B. 

ii. Compute a new basis set (j)i{T) by using the relation 0(r) = M~^^(t), then compute 
b = [61, 62, ... , b^y by bi = (j)i{T)dT. 

4. Let the matrix A'^ = [aij] contain the coefRcients of the SPRK tableau. As before, there 
are two way to obtain A'^: 

i. Compute coefficients A^ = [Aij], where Aij = J^"^ ipj{T)dT. Then, the matrix is 
given by A"^ = [o^] = A'f'M''^. 

ii. Compute A'^ = [aij], where a,j = Jq' (l)j{T)dT, directly by using the new basis 
functions </>(•) = M-V(-)- 

5. Compute the coefRcients A'^ = [dij\ by using a^- = bh^i^^ 



4.5. Examples. In this subsection, we will consider four examples to illustrate the above procedure for 

constructing variational integrators. 



Example 2. We consider one-stage methods. Choose the basis function tpi = 1; then for any quadrature 
point c\, the matrix M = [V'i(ci)] = [1] is invertible. 

i. If ci = 0, then bi = l,au = 0,aii = 1, which is the symplectic Euler method. 

ii. If ci = ^, then bi = l,aii = ^,aii = i, which is the midpoint rule. 

iii. If ci = 1, then bi = I, an ~ I, an = 0, which is the adjoint symplectic Euler method. 



Example 3. Choose basis functions tpi = l,tp2 = cos{ttt). If we choose quadrature points Ci = 0,C2 = 1, 
then we obtain 



M 



One can easily compute 



V'l(ci) V'l(C2) 
V'2(ci) V2(C2) 





1 1 " 




1 -1 



and = - 

2 



1 1 
1 -1 



Bi= [ MT)dT = l, B2= [ 

Jo Jo 



Therefore we get b = [61, 62]"^ = M ^B = [^, ^]'^ , which is the trapezoidal rule. We also compute 



A'/' 



J^'Mr)dr CMr)d 








0" 




1 






Therefore, the matrix 



lo' MT)dT /o>2(r)dT 
j;;'Mr)dr J^'Mr)dT^ 



A^M-^ = 





1 1 

2 2. 



By using the relationship cii 



— , / ^ , one obtains 





1 1 
1 2. 



Thus, we obtain the Stdrmer-Verlet method. Interestingly, the Stdrmer-Verlet method is typically derived as 
a variational integrator by using linear interpolation, i.e., ipi = 1,V'2 = t, and the trapezoidal rule. 



Example 4. // we choose basis functions ipi = l,ijj2 = cos(7rT),^3 = sin{'KT) and quadrature points c\ 



0,C2 



,C3 



1, we obtain a new method which is second-order accurate, and the coefficients of the SPRK 



method are given by 
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1 j_ 

4 TT 

ir-2 2 



4ir 
7r-2 



ir-2 



7r-2 



7r^-2ir-4 


7r-2 


TT — 4 





27r^-47r 


2ir 


ir^-2ir 


1 


7r-2 


j_ 





2 


2ir 


TT 


7r^-27r+4 


7r-2 


1 





27r2-47r 


27r 


7r-2 




7r-2 


2 






27r 


7r 


2-K 



Example 5. Chebyshev quadrature (see p. 415 of '14' j is designed to approximate weighted integrals of the 
form 



f f{x)w{x)dx = bY,f{^i) + E[f{x)l 
■^^1 i=i 



with an equally weighted sum of the function values at the quadrature points Xi, and an error term E[f{x)\. 
The weight b is chosen so that the quadrature is exact for f{x) = 1, i.e., b — ^ J^-^ w{x)dx. We are primarily 
interested in the case when the weight function w{x) ~ 1, in which case the quadrature formula becomes 



where the quadrature points Xi are given by the roots of polynomials (see p. ^18 of [TT j, the first three of 
which are given by 



(44) 



Go(a;) = l, Gi(x) 



G2(a:) = i(3x2 - 1), G3(x) - i(2a;3-x). 



The error term associated with the s -point formula is given by 



E 



(s+l)! 
(s+2)! 



s odd, 
8 even. 



where 



J^-^xGs{x)dx s odd, 
J^-^ x^Gs{x)dx s even. 



The error term implies that the quadrature has degree of precision s for odd s and degree of precision s + 1 
for even s. Note that the roots Xi of the polynomials Gi are in the interval [—1,1], so after a change 
of coordinates, we obtain quadrature points Ci in the interval [0,1]. Then, we use Lagrange polynomials 
associated with these quadrature points to construct variational integrators for s — 1,2, 3. 

i. one-stage, second- order method 



1 1 

2 2 



ii. two-stage, fourth-order method 



1 V3 

2 6 

1 -i_ VI 

2 6 



1 

4 

1 J_ vi 

4 6 



V3 
6 



iii. three-stage, fourth-order method 

V2 



1 _ \/2 

2 4 
1 

2 

1 _L v2 

2 4 



1 

6 ' 48 

1 -|_ %/I 

6 8 

1 I 5\/2 
6 48 6 



1 _ V2 1 _ 5\/2 

6 6 6 48 
1 1 V2 

6 6 8 

1 j_ VZ 1 _ 5^/2 

6 6 48 



1 


1 






2 


2 








1 






1 


V3 


1 


1 V3 


2 


6 


4 


4 6 


1 


1 


1 1 V3 


1 


2 


6 


4 ' 6 


4 






1 


1 






2 


2 



1 _ V2 

2 4 
1 

2 

1 _L VI 

2 4 



s/2 1 _ \/2 1 _ 5^2 

48 6 8 6 48 

1 1 V2 

6 6 6 

1 _L 5\/2 1 I \/2 1 I V2 

6"*" 48 6~*'8 6"^ 48 



1 -|_ vi 

6 6 



For s = 1,2, we obtain the same methods as the ones obtained using Gauss-Legendre quadrature, which 
are the midpoint rule, and the two-stage, fourth-order method, respectively. For s = 3, we obtain a three- 
stage SPRK that is fourth-order. The order of the SPRK methods above were determined using partitioned 
Runge-Kutta order theory [5] . 
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5. Momentum Preservation and Invariance of the Discrete Right Hamiltonian Function 
Momentum Maps. First, we recall the definition of a momentuni map defined on T*Q given in [T]. 

Definition 1. Let {PtLo) he a connected symplectic manifold and ^ : G x P ^ P be a symplectic action 
of the Lie group G on P, i.e., for each g G, the map <^g : P P; x ^ ^{g^x) is symplectic. We say 
that a map J : P ^ q* , where g* is a dual space of the Lie algebra q of G, is a momentum map for the 
action <f> if for every ^ G 0, dJ{S^) = i^pUJ, where J(i^) : P — ^ K is defined by J{S.){x) — J{x) ■ and is 
the infinitesimal generator of the action corresponding to ^. In other words, J is a momentum map provided 
^J(5) = all ^ eg. 

For om purposes, we are interested in the case P = T*Q, and uj — dq^ A dpi is the canonical symplectic 
two-form on T*Q. This gives a momentum map of the form J : T*Q — > g*, and we describe the construction 
given in Theorem 4.2.f0 of 1 . Notice that lu is exact, since uj — --dO ~ —d{pidq^). Consider an action 
$g that leaves the Lagrange one-form 9 invariant, i.e., $*6' = 6 for all g E G. Then, the momentum map 
J : T*Q 0* is given by 

(45) J{x)-^^t^p0{x). 

We can show that this satisfies the definition of the momentum map given above by using the fact that $g 
leaves the one-form 9 invariant for all g G G, and is the infinitesimal generator of the action corresponding 
to ^. This implies that the Lie derivative of along the vector field ^p vanishes, i.e., £^p9 — 0, for all ^ G g. 
By Cartan's magic formula, 

= £^p9 ~ i^pd9 + di^p9, 

therefore, di(^p9 — —i^pd9 — i^pLo. As such, J{£^){x) = i^p9 satisfies the defining property, dJ{£,) = i^pOJ, of 
a momentum map. Then, J{x) ■ ^ = J{£,){x) = i^p9{x). Moreover, by Theorem 4.2.10 of [1], this momentum 
map is Ad*-equivariant. 

Let $ : G X Q — >■ Q be an action of the Lie group G on Q. We will give the coordinate expression for the 
cotangent lifted action $^ . In coordinates, we denote $g-i : Q ^ Q by = <I>^_i((5), then its cotangent 
lifted action <^'^'Q : G x T*Q ~> T*Q is given by 

(46) ^^'Q{g,q,p)=T*<i>g-.{q,p) ^ {^^^{q),p^-^^ , 

where T*$g-i means the cotangent lift of the action $3-1. 

In the following proposition, we give the coordinate expression for the cotangent lifted action, and show 
that it leaves the Lagrange one- form 9 — pid(f invariant. 

Proposition 3. Given an action ^ : G x Q of a Lie group G on Q, the cotangent lifted action ; 
G X T*Q — > T*Q leaves the Lagrange one-form 9 = ptdq' invariant. 

Proof. Given g E G, let the cotangent lifted action of g on {q,p) be denoted by {Q,P) — $J '^{q,p), the 
components of which are given by = $g(q) and Pi — Pj-Q§Tj;^- Then, a direct computation yields 

da^ d^Uq) 

(47) p^dQ^ =P^d<I>;(,) ^p^dq^. 

This shows that <I>J leaves the Lagrange one- form pidq^ invariant. □ 

Corresponding to the cotangent lift action ^, for every ^ G 0, the momentum map J : T*Q ^ q* 
defined in (j45p has the following explicit expression in coordinates. 



(48) Jiaq) ■ C = i^p{pidq'){aq) ^p- ^giq) = P ' ^ 

where a, — {q,p) G T*Q. 



*cxp(.c)(g'), 

e=0 
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Discrete Noether's Theorem. In the discrete case, consider the one-step discrete flow map i^^-t 
{qo,Po) I— > {qi,pi) defined by the discrete right Hamihon's equations, 

(49) po = DiH+{qo,p,), q,^D2H+{qo,pi). 



We win show that if the generalized discrete Lagrangian, Rd{qo, qi,Pi) = Piqi—H^{qo,pi), is invariant under 
the cotangent lifted action, then we have discrete momentum preservation, which is the discrete analogue of 
Noether's theorem for discrete Hamiltonian mechanics. 

Theorem 4. Let ^ be the cotangent lift action of the action $ on the configuration manifold Q. If 
the generalized discrete Lagrangian Rdiqo, qi,Pi) ~ Piqi — H'^ {qo , pi) is invariant under the cotangent lifted 
action ^, then the discrete flow map of the discrete right Hamilton's equations preserves the momentum 
map, i.e., F*^J— J. 



Proof. In coordinates, let ((7o;-Po) 
Piqi — {qo,pi), we have that 



Cp?.o('?o,Po) and (qlpl) ^p^o 



(91 J Pi)- From the invariance of 



(50) 



d 

'-de 



€=0 





d 


= P1 


'de 




d 


= P1 






d 


= P1 


d'e 



Ht{qlp\)] 
d 



e=0 



p\-D,H2{qo,Pi) -r 



e=0 



9i +91 



e=0 



de 



Pi -Po 



e=0 



e=0 

Pi •Cq('7i) 



*cxp(eC)(gi) - Po 

Po •Cq('7o), 



de 



de 



de 

=0 

oxp(eC)(9o) 



Qo 



e=0 

d_ 
d'e 



D2H+{qo,pi) - 



Pi 



c=0 



Pi 



e=0 



e=0 



where we used the discrete right Hamilton's equations (|49l) in going from the second to the third line, and 
^o = ^oxp(e^)('7o) and qf = 'I*cxp(e^) (91) are used in going from the third to the fourth line. Then, by the 
definition of F„+ and J, ([50|) states that F* + J — J. 



□ 



G-invariant generalized discrete Lagrangians from G-equivariant interpolants. We now provide 
a systematic means of constructing a discrete Hamiltonian, so that the generalized discrete Lagrangian 
Rdiqo,qi,Pi) = Piqi — H^{qo,pi) is G-invariant, provided that the generalized Lagrangian R{q,q,p) = 
pq — H{q,p) is G-invariant. 

Our construction will be based on an interpolatory function ip : — ?> G^([0, h],Q), that is parameterized 
by r-f 1 internal points q'^ G Q, defined at the times = doh < dih < ■ ■ ■ < drh < h, i.e., ip {dj^h; {g'^jy^o) = 
q^'. We also use a numerical quadrature formula given by quadrature weights bi and quadrature points c^. 
We denote the momentum at the time Cih by p^ . Then, we construct the following discrete Hamiltonian, 



(51) 



Hai<lo,Pi) 



ext 

9" =90 



PI • ^ (/i; {q''}l=o) - E (^^ (^''^^ W}1=^) y 



where R{q, q,p) = pq ^ H{q,p). An interpolatory function is G-equivariant if 

V{t-A9q'^]l=,)=gv{t-Aq''}U,). 

Then, a G-invariant discrete Hamiltonian can be obtained if we use G-equivariant interpolatory functions. 

Lemma 3. Let G be a Lie group acting on Q, such that gQ—Q, for all g G G. If the interpolatory function 
y^i^j {9q'^}l^=o) ^'^ G-equivariant, and the generalized Lagrangian R : TQ © T*Q — > M, 

R{q,q,p) =pq- H{q,p), 

is G-invariant, then the generalized discrete Lagrangian Rd : Q x T*Q ^ M., given by 

Rd{qo,qi,Pi) =Pigi - Hj{qo,pi), 
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where, 



Hdiqo,Pi) = ext 



Pi ■ V (h; {qT.=o) ~ E (^V' ^^^-^ WM=^) ,P\ 



is G -invariant. 



Proof. To streamline the notation, we denote the cotangent hfted action of G on Q by $J '^{q,p) — {gq, gp). 
First, we note that 



Rd{qQ,qi,pi) = piqi - ext 

9" =90 



ext 

9" =90 



Pi ■ V {h- {q''}^^) - E hR {Tif ich; {q'-r^^o) ,f 

i— 

Y,hR{T^{c,h-{q''}l^^),p\ 



Then, 



Rd{gqo,gqi,gpi) ^ ext 

9° =990 



ext 

91° =910 



ext 

9"eQ,p'ec 

9^=90 

ext 

9"eQ,p'ec 

9^=90 



J2hR{T^ic,h;{r}:=o),p') 

S 

J2hR{Tv{c,h;{gq-'}:^o).9P') 

i=l 

J2b,R{TLg-T^ic.h;{qTu^,),gp^) 

s 

Y,hR{T^{c,h-{q'^}l^,),p') 



= Rd{qo,qi,Pi), 

where we used the identification q'^ = gq'^ in the second equahty, the G-equivariance of the interpolatory 
fimction and the property that gQ = Q in the third equahty, and the G-invariance of the generahzed 
Lagrangian in the fourth equahty. □ 

In view of Theorem 01 and the above lemma, if we use a G-equivariant interpolatory function to construct 
a discrete Hamiltonian as given in (|5ip. then the discrete flow given by the discrete right Hamilton's equa- 
tions will preserve the momentum map J : T*Q — s- g*. 



Natural Charts and G-equivariant interpolants. Following the construction in [T5], we use the group 
exponential map at the identity, expg : g — i> G, to construct a G-equivariant interpolatory function, and a 
higher-order discrete Lagrangian. As shown in Lemma [31 this construction yields a G- invariant generalized 
discrete Lagrangian if the generalized Lagrangian itself is G- invariant. 

In a finite-dimensional Lie group G, exp^ is a local diffeomorphism, and thus there is an open neighborhood 
[/ C G of e such that exp~^ : U u C q. When the group acts on the left, we obtain a chart ipg : Lgll u 
at 5 e G by 

'ipg = exp^^ °Lg-i. 
Lemma 4. The interpolatory function given by 

is G-equivariant. 
Proof. 

ifigg^lTh) = ^^glQ',[j2^=o^99'>i99''%,'^iT)j 
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= L,gO exp, (^^^^exp-i((5g0)-i(5g''))/;,,(r)) 

= Lgifiig^iTh). □ 

This G-equivariant interpolatory function based on natural charts allows one to construct discrete Lie 
group Hamiltonian variational integrators that preserve the momentum map. 

6. Conclusions and Future Directions 

In this paper, we provided a variational characterization of the Type II generating function that generates 
the exact flow of Hamilton's equations, and show how this is a Type II analogue of Jacobi's solution of 
the Hamilton-Jacobi equation. This corresponds to the exact discrete Hamiltonian for discrete Hamiltonian 
mechanics, and Galerkin approximations of this lead to computable discrete Hamiltonians. In addition, we 
introduced a discrete Type II Hamilton-Jacobi equation, which can be viewed as a composition theorem for 
discrete Hamiltonians. 

We introduced generalized Galerkin variational integrators from both the Hamiltonian and Lagrangian 
approach, and when the Hamiltonian is hyperregular, these two approach are equivalent. Furthermore, we 
demonstrated how these methods can be implemented as symplectic partitioned Runge-Kutta methods, 
and derived several examples using this framework. Finally we characterized the invariance properties of a 
discrete Hamiltonian which ensure that the discrete Hamiltonian flow preserves the momentum map. 

We are interested in the following topics for future work: 

• Lie-Poisson Reduction and Connections to the Hamilton-Pontryagin principle. Since we provided 
a method for constructing discrete Hamiltonians that yields a numerical method that is momentum 
preserving, it is natural to consider discrete analogues of Lie-Poisson reduction. In particular, the 
constrained variational formulation of continuous Lie-Poisson reduction [6] appears to be related 
to the Hamilton-Pontryagin variational principle 27 . It would be interesting to develop discrete 
Lie-Poisson reduction [18] from the Hamiltonian perspective, in the context of the discrete Hamilton- 
Pontryagin principle [16| I25j . 

• Extensions to Multisymplectic Hamiltonian PDEs. Multisymplectic integrators have been developed 
in the setting of Lagrangian variational integrators |17j , and Hamiltonian multisymplectic integra- 
tors [4]. In the paper [20], the Lagrangian formulation of multisymplectic field theory is related 
to Hamiltonian multisymplectic field theory [3]. It would be interesting to construct Hamiltonian 
variational integrators for multisymplectic PDEs by generalizing the variational characterization of 
discrete Hamiltonian mechanics, and the generalized Galerkin construction for computable discrete 
Hamiltonians, to the setting of Hamiltonian multisymplectic field theories. 
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